--- title: Graphics keywords: fastai sidebar: home_sidebar summary: "This module encapsulates the functions used for displaying an optimisation function" description: "This module encapsulates the functions used for displaying an optimisation function" nb_path: "01_graphics.ipynb" ---
{% raw %}
{% endraw %} {% raw %}

contour[source]

contour(function:Ifunction, x_min=-5, x_max=5, y_min=-5, y_max=5, mesh_size=100)

Returns a (x, y, z) 3D coordinates, where z = function(x,y) evaluated on a mesh of size (mesh_size, mesh_size) generated from the linear space defined by the boundaries returned by function.domain().

This function is usually used for displaying the contour of the diven function.

{% endraw %} {% raw %}
{% endraw %} {% raw %}

zoom[source]

zoom(x_domain, y_domain, zoom_factor)

{% endraw %} {% raw %}
{% endraw %}

The idea of zooming is actualy closely related to the min and max boundaries set on the plots:

  • if we reduce the distance between them, we zoom in
  • if we increase the distance, we zoom out

Since zoom is usually a proportion of the current view, there's a bit of math to do in order to get the computations correct.

This section is highly specific and beside the main point but I like to keep it just because it's time spent that I don't want to be lost to history. If you're reading this, you might probbably want to skip it since it isn't clever nor informative..

So the main idea of the zooming that I'm going to implement is described by the image bellow:

  • we have some set 'min' and 'max' values
  • we have a mean (center) between them
  • we want to move both min and max closer to the mean by the same amunt (percentage), f called a zooming factor.

image.png

{% raw %} 
assert zoom((-5, 5), (-5, 5), 0.1) == ((-4.5, 4.5), (-4.5, 4.5))
{% endraw %}

Funny enough, these can be vectorized by the numpy transformations bellow

{% raw %} 
d = np.array([(-5, 5), (-5, 5)])
zoom_factor = 0.1

means = d.mean(axis=-1)             # computing the means
distances = np.abs(d - means)
change = distances * zoom_factor    # compute the change need 
change_with_direction = change * [1, -1] # add signs for the direction of changes (mins should increase, maxes should decrese in value)
zoomed_d = d + change_with_direction

d + np.abs(d - d.mean(axis=-1)) * [1, -1] * zoom_factor # single line transformation

array([[-4.5,  4.5],
       [-4.5,  4.5]])
{% endraw %}

Now, we see that the above computation does, an abs where we eliminate the sign, and right after that we add it back by multipligin with [-1, 1]. If we think in terms of moving from the mean left and right we can simplify the formula a bit, as follows:

{% raw %} 
d.mean(axis=-1) - (d.mean(axis=-1) - d) * (1 - zoom_factor)

array([[-4.5,  4.5],
       [-4.5,  4.5]])
{% endraw %}

If we analitically decompose the (1 - 0.1) term and simplify the result, we remain with the bellow formula:

{% raw %} 
d + (d.mean(axis=-1) - d) * zoom_factor

array([[-4.5,  4.5],
       [-4.5,  4.5]])
{% endraw %}

Which is almost identical with the one we've started from but, as we've observed, does not trim the signs and adds them back.

3D Plot

{% raw %}

__log_contour_levels[source]

__log_contour_levels(zz, number_of_contour_lines=35)

{% endraw %} {% raw %}
{% endraw %} {% raw %}

plot_function_3d[source]

plot_function_3d(function:Ifunction, ax=None, azimuth=45, angle=45, zoom_factor=0, show_projections=False, contour_log_scale=True)

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
from optimisations.functions import himmelblau
fig = plt.figure(figsize=(13, 10))
ax = plt.axes(projection='3d')
plot_function_3d(himmelblau(), ax=ax, azimuth=20, angle=225)
{% endraw %}

Let's make it interactive so we can play with it a bit.

Note: This code will probably not work in the blog post

{% raw %} 
import ipywidgets as widgets
from ipywidgets import interact, interact_manual

@interact
def plot_interactive(azimuth=45, rotation=45, zoom_factor=(-1, 1, 0.1), show_projections={True, False}):
    return plot_function_3d(function=himmelblau(), azimuth=azimuth, angle=rotation, zoom_factor=zoom_factor, show_projections=show_projections)
{% endraw %}

2D plot

Contour with angle rotation options

Note that we have a 3D plot what can handle rotations, we need to allow this capability to the 2D plot as well, since we wish the two charts to move in sync.

What we need to do is use linear algebra and rotate the initial meshgrid of points. We accomplish this by simply multiplying it with a rotation matrix.

Adding rotation to the second plot

{% raw %} 
angle = 40
xx, yy, zz = contour(himmelblau())


radians = angle * np.pi/180

# counter-clockwise rotation matrix
# https://stackoverflow.com/questions/29708840/rotate-meshgrid-with-numpy
rotation_matrix = np.array([
    [np.cos(radians), -np.sin(radians)],
    [np.sin(radians), np.cos(radians)]
])
xx, yy = np.einsum('ji, mni -> jmn', rotation_matrix, np.dstack([xx, yy]))

plt.contourf(xx, yy, zz, levels=1000, cmap='Spectral', norm=colors.Normalize(vmin=zz.min(), vmax=zz.max()), alpha=0.3)

<matplotlib.contour.QuadContourSet at 0x13706fdd8>
{% endraw %}

We've successfully rotated the contour, but we see that the axis were left unchanges. So we didn't rotate the plot, but merely it's contents.

At this point, I guess there are three options:

  • use some image processing packages (like PIL, imagemagick), export the 2D plot as an image, rotate it and then display it
  • hide the axes so we don't see that we've actually rotated only the content.
  • use a 3D plot (which can easily suport rotation as we've seen) and set a perpendicular viewing angle (birds-eye view) so the plot looks like a 2D one.

The first option looks like the biggest hack of all since it involves adding at least 2 new dependencies for this sole purpose (plus the additional computation).

The last one is easy to see without further experiments, so we only need to see how option 2 looks like.

Drawing the 2D plot on a 3D Axis and viewed from above. This makes it possible to dispay the rotated axis as well but also makes the plot smaller.

{% raw %} 
fig = plt.figure()

ax_ = fig.add_subplot(1, 2, 2, projection='3d')

angle = 125

## plot_function_2d
xx, yy = np.meshgrid( 
    np.linspace(-5, 5, num=100),
    np.linspace(-5, 5, num=100)
)
zz = himmelblau()(xx, yy)

ax_.contour(xx, yy, zz, levels=200, offset=0, cmap='Spectral', norm=colors.Normalize(vmin=zz.min(), vmax=zz.max()), alpha=0.5, zorder=1)
ax_.view_init(90, angle)

plt.tight_layout()

/Users/cristi/Envs/optimisations/lib/python3.6/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)
{% endraw %}

I think I will stick with the 2D approach as the contour is bigger.

Rotation function

The rotate function bellow looks way more complex than we've sketched it out above, and this is because it is now handling two use cases that it needs to disambiguate:

  • the case when it receives meshgrids (xx, yy) when called to rotate the contour plots & the case when it receives points we want to draw on the contour (like the minimum values or the trace of an optimizer) which are not of the same shape as a meshgrid so have a diffrent type of math operation applyed on them (mainly due to vectorisation).
{% raw %}

rotate[source]

rotate(_x:ndarray, _y:ndarray, angle=45)

{% endraw %} {% raw %}
{% endraw %}

The np.einsum part is the core of the function, and it's operation is more complex, leaving its explanation for a future post.

Another interesting bit of it is the way the code bellow works, namely how the parameters are returned.

xx, yy = np.einsum('ij, mnj -> imn', rotation_matrix, np.dstack([xx, yy]))

Normally, any np. prefixed operation returns a np.ndarray but in this case we see that the enisum function is able to destructure the return into two separate parameters xx and yy. This happens with the help of the np.dstack function.

Let's see what it does, bellow:

{% raw %} 
print(np.dstack([xx, yy]).shape)
np.dstack([xx, yy])[:3, :3, :]

(100, 100, 2)

array([[[-5.       , -5.       ],
        [-4.8989899, -5.       ],
        [-4.7979798, -5.       ]],

       [[-5.       , -4.8989899],
        [-4.8989899, -4.8989899],
        [-4.7979798, -4.8989899]],

       [[-5.       , -4.7979798],
        [-4.8989899, -4.7979798],
        [-4.7979798, -4.7979798]]])
{% endraw %}

So dstack just made a stack of it's arguments, on a new, rightmost axis (the last 2 value from the shape). This shape, coupled with the specified einsum operation, where both i and j are equal to 2 gives a result with the shape of (2, 100, 100) semantically similar to a tuple (xx.shape(100, 100), yy.shape(100, 100)). This shape enables the python language to destructure this return into 2 independent values, the ones we actually wanted.

Single points rotation

The gist of single points rotation is a matrix multiplication but since the meshgrid is implemented in einsum notation as is rather elegant, at least because of the arguments decomposition trick, I'll try experimenting a bit to replicate the matrix multiplication it with einsum as well.

{% raw %} 
angle = 225

# apply rotation
angle = (angle) % 360
radians = angle * np.pi/180


# clockwise rotation matrix
rotation_matrix = np.array([
    [np.cos(radians), -np.sin(radians)],
    [np.sin(radians), np.cos(radians)]
])


points = np.array([
    [1, 1],
    [-2, 2],
    [3, -3],
    [4, 4]
])

function = himmelblau()
coords = function.coord(points)
coords[:, [0, 1]] @ rotation_matrix

array([[-1.41421356e+00, -2.22044605e-16],
       [ 4.44089210e-16, -2.82842712e+00],
       [-4.44089210e-16,  4.24264069e+00],
       [-5.65685425e+00, -8.88178420e-16]])
{% endraw %}

Ok, this is the (correct) result we're getting using the plain matrix multiplication. Our goal is the einsum operation that computes these, but reshapes them as well in (2, 4) so we can later on decompose them as (4,) and (4,) into x and y coordinates.

We start of with 4 points, where we also have the z value (the function evaluation on the first two coordinates), so a shape of (4, 2+1)

{% raw %} 
coords.shape

(4, 3)
{% endraw %}

We'll be receiving x and y values in the following shape, so this is what we have to work with:

{% raw %} 
_x = coords[:, [0]]
_y = coords[:, [1]]
_x.shape, _y.shape

((4, 1), (4, 1))
{% endraw %}

Since the rotation_matrix is of shape (2,2), the final einsum operation is:

mi, ij -> jm

or 
(nr_points, 2), (2, 2) -> (2, nr_points)
{% raw %} 
__x, __y = np.einsum('mi, ij -> jm',  np.hstack((coords[:, [0]], coords[:, [1]])), rotation_matrix)
__x, __y

(array([-1.41421356e+00,  4.44089210e-16, -4.44089210e-16, -5.65685425e+00]),
 array([-2.22044605e-16, -2.82842712e+00,  4.24264069e+00, -8.88178420e-16]))
{% endraw %}

Actual function

Again, let's wrap everything up in a single function to see where we are up until now.

{% raw %}

plot_function_2d[source]

plot_function_2d(function:Ifunction, ax=None, angle=45, zoom_factor=0, contour_log_scale=True)

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
plot_function_2d(himmelblau(), contour_log_scale=False)
{% endraw %}

Joint plot

Using both 2D and 3D plotting function we can combine them to be shown on the same figure, and using the same rotation.

As allways, we will first sketch the code in bulk then tie it up in a single function.

{% raw %} 
fig = plt.figure(figsize=(26, 10))

ax_3d = fig.add_subplot(1, 2, 1, projection='3d')
ax_2d = fig.add_subplot(1, 2, 2)

angle = 225

function = himmelblau()

plot_function_3d(function, ax=ax_3d, azimuth=30, angle=angle)
plot_function_2d(function, ax=ax_2d, angle=angle)
{% endraw %}

The 2D and 3D plots can now be combined in a unified function that you can see bellow:

{% raw %}

plot_function[source]

plot_function(function:Ifunction, angle=45, zoom_factor=0, azimuth_3d=30, fig=None, ax_2d=None, ax_3d=None, contour_log_scale=True)

{% endraw %} {% raw %}
{% endraw %}

Usage examples

{% raw %} 
from optimisations.functions import *
plot_function(himmelblau(), angle=225, contour_log_scale=False)

(<Figure size 1872x720 with 2 Axes>,
 <matplotlib.axes._subplots.Axes3DSubplot at 0x138039dd8>,
 <matplotlib.axes._subplots.AxesSubplot at 0x138068470>)
{% endraw %} {% raw %} 
plot_function(mc_cormick(), angle=225, contour_log_scale=False)

WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

(<Figure size 1872x720 with 2 Axes>,
 <matplotlib.axes._subplots.Axes3DSubplot at 0x1385665c0>,
 <matplotlib.axes._subplots.AxesSubplot at 0x137579c18>)
{% endraw %} {% raw %} 
plot_function(holder_table(), angle=225, contour_log_scale=False)

(<Figure size 1872x720 with 2 Axes>,
 <matplotlib.axes._subplots.Axes3DSubplot at 0x1362c1f60>,
 <matplotlib.axes._subplots.AxesSubplot at 0x13a6f6f60>)
{% endraw %}

Let's see all of our functions defined

{% raw %}

plot_all_functions[source]

plot_all_functions(functions:dict)

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
from optimisations.functions import Function
plot_all_functions(Function)

/Users/cristi/Envs/optimisations/lib/python3.6/site-packages/ipykernel_launcher.py:13: UserWarning: Log scale: values of z <= 0 have been masked
  del sys.path[0]
/Users/cristi/Envs/optimisations/lib/python3.6/site-packages/ipykernel_launcher.py:13: UserWarning: No contour levels were found within the data range.
  del sys.path[0]
{% endraw %}